ReactionCO2SYS

Carbonate chemistry using PALEOcarbchem.

Two different pH solver approaches are supported, set by Parameter solve_pH:

• solve: iterative solution for pH using an internal Newton solver, at given TAlk (eg TAlk is a StateExplicit Variable, using explicit ODE solver)
• speciation: calculate speciation given pHfree.
• speciationTAlk: calculate speciation given pHfree, also add alkalinity contributions to TAlk_calc (eg for use in combination with a ReactionConstraintReservoir that provides a total variable TAlk and H (as pHfree) as a primary species, for use with a DAE solver to solve the algebraic constraint for pHfree at a given TAlk).

NB: Two options have been removed from Parameter solve_pH:

• constraint: additional State Variable for pH and algebraic Constraint for TAlk (requires DAE solver). [replace with speciationTAlk and a ReactionConstraintReservoir]
• implicit: TAlk is a Total Variable, a function of a State Variable for pH (requires DAE solver) [replace with speciationTAlk and a ReactionImplicitReservoir]

Parameters

• WhichKs[Int64]=10, default_value=10, description="PALEOcarbchem WhichKs choice of eqb constant data"
• WhoseKSO4[Int64]=1, default_value=1, description="PALEOcarbchem WhoseKSO4 choice of eqb constant data"
• components[Vector{String}]=["Ci", "S", "F", "B", "Omega"], default_value=["Ci", "S", "F", "B", "Omega"], allowed_values=["Ci", "S", "F", "P", "B", "Si", "H2S", "NH3", "Omega"], description="PALEOcarbchem optional contributions to TAlk"
• defaultconcs[Vector{String}]=["TS", "TF", "TB", "Ca"], default_value=["TS", "TF", "TB", "Ca"], allowed_values=["TS", "TF", "TB", "Ca"], description="PALEOcarbchem default concentrations (modern values calculated from salinity)"
• outputs[Vector{String}]=["pCO2", "xCO2dryinp"], default_value=["pCO2", "xCO2dryinp"], allowed_values=["H", "OH", "TS", "HSO4", "TF", "HF", "TCi", "CO2", "HCO3", "CO3", "CAlk", "fCO2", "pCO2", "xCO2dryinp", "TB", "BAlk", "TP", "H3PO4", "H2PO4", "HPO4", "PO4", "PAlk", "TSi", "SiAlk", "TH2S", "H2S", "HSAlk", "TNH3", "NH4", "NH3Alk", "Ca", "OmegaCA", "OmegaAR"], description="PALEOcarbchem output concentrations etc to include as output variables"
• output_pHtot[Bool]=true, default_value=true, description="true to output pHtot (requires TS_conc and TF_conc)"
• solve_pH[String]="solve", default_value="solve", allowed_values=["solve", "speciation", "speciationTAlk"], description="'solve' to solve iteratively for pH, 'speciation' to just provide speciation at supplied [pHfree], 'speciationTAlk' to provide speciation and add to TAlk_calc"
• pHtol[Float64]=2.220446049250313e-14, default_value=2.220446049250313e-14, description="with parameter solve_pH=solve only (internal Newton solver): pH tolerance for iterative solution"
• simd_width[String]="1", default_value="1", allowed_values=["1", "FP64P2", "FP64P4", "FP32P4", "FP32P8"], description="with parameter solve_pH=solve only (internal Newton solver): use SIMD ("1" - no SIMD, "FP64P4" - 4 x Float64, etc)"

Methods and Variables for default Parameters

• do_modern_default_concs
  • volume (m^3), VT_ReactDependency, description="Cell volume"
  • rho_ref (kg m^-3), VT_ReactDependency, description="density conversion factor"
  • temp (K), VT_ReactDependency, description="temperature"
  • pressure (dbar), VT_ReactDependency, description="pressure"
  • sal (psu), VT_ReactDependency, description="salinity"
  • TS_conc (mol m^-3), VT_ReactProperty, description="TS_conc modern default scaled from salinity"
  • TF_conc (mol m^-3), VT_ReactProperty, description="TF_conc modern default scaled from salinity"
  • TB_conc (mol m^-3), VT_ReactProperty, description="TB_conc modern default scaled from salinity"
  • Ca_conc (mol m^-3), VT_ReactProperty, description="Ca_conc modern default scaled from salinity"
• do_carbchem
  • volume (m^3), VT_ReactDependency, description="Cell volume"
  • rho_ref (kg m^-3), VT_ReactDependency, description="density conversion factor"
  • temp (K), VT_ReactDependency, description="temperature"
  • pressure (dbar), VT_ReactDependency, description="pressure"
  • sal (psu), VT_ReactDependency, description="salinity"
  • pHfree (), VT_ReactProperty, description="-log10(hydrogen ion concentration), this is used to keep previous value as a starting value for internal Newton solver"
  • TAlk_conc (mol m^-3), VT_ReactDependency, description="total TAlk concentration"
  • pHtot (), VT_ReactProperty, description="-log10([H] + [HS]), pH on total scale"
  • TCi_conc (mol m^-3), VT_ReactDependency, description="TCi_concinput concentration"
  • TS_conc (mol m^-3), VT_ReactDependency, description="TS_concinput concentration"
  • TF_conc (mol m^-3), VT_ReactDependency, description="TF_concinput concentration"
  • TB_conc (mol m^-3), VT_ReactDependency, description="TB_concinput concentration"
  • Ca_conc (mol m^-3), VT_ReactDependency, description="Ca_concinput concentration"
  • pCO2 –> %reaction%pCO2 (atm), VT_ReactProperty, description="CO2 partial pressure (fugacity corrected)"
  • xCO2dryinp –> %reaction%xCO2dryinp (), VT_ReactProperty, description="mixing ratio of CO2 in dry air at 1 atm (always > pCO2 due to H2O vapour pressure)"

PALEOcarbchem

Carbonate chemistry equilibrium translated from Matlab CO2SYS v1.1 (van Heuven et al., 2011), (Lewis and Wallace, 1998) and refactored for speed.

Implementation based on CO2SYS constants, with reimplemented equilibrium calculation.

NB:

• Although effective equilibrium constants can be calculated for different pH scales, TAlk or pH calculations using PALEOcarbchem.calculateTAfromTCpHfree!

require constants to be on the free pH scale.

• Effective equilibrium constants are empirical measurements hence are only accurate for modern seawater composition.

Example usage with default choice of constants, all components enabled:

julia> PALEOcarbchem.ComponentsAllStrings  # all available components
("Ci", "S", "F", "B", "P", "Si", "H2S", "NH3", "Omega")

julia> comps, concinputs = PALEOcarbchem.get_components_inputs(["Ci", "S", "F", "B", "P", "Si", "H2S", "NH3", "Omega"])
((Ci = Val{true}(), S = Val{true}(), F = Val{true}(), B = Val{true}(), P = Val{true}(), Si = Val{true}(), H2S = Val{true}(), NH3 = Val{true}(), Omega = Val{true}()), ["TCi", "TS", "TF", "TB", "TP", "TSi", "TH2S", "TNH3", "Ca"])

julia> println(concinputs) # concentrations required (in addition to TAlk)
["TCi", "TS", "TF", "TB", "TP", "TSi", "TH2S", "TNH3", "Ca"]  

julia> options = (; WhichKs=Val(10), WhoseKSO4=Val(1), pHScale=Val(3), Components=comps)
(WhichKs = Val{10}(), WhoseKSO4 = Val{1}(), pHScale = Val{3}(), Components = (Ci = Val{true}(), S = Val{true}(), F = Val{true}(), B = Val{true}(), P = Val{true}(), Si = Val{true}(), H2S = Val{true}(), NH3 = Val{true}(), Omega = Val{true}()))

julia> C = zeros(length(PALEOcarbchem.ConstNames));

julia> PALEOcarbchem.calc_constants!(C, 25.0, 1000.0, 35.0, Options=options) # Temp(C), P(dbar), Sal(psu)

julia> C_NT = NamedTuple{PALEOcarbchem.ConstNames}(C);  # convert to NamedTuple for convenience

julia> map(x -> @sprintf("%.14e", x), C_NT)
(TF = "6.83258396883673e-05", TS = "2.82354341328601e-02", fH = "7.13404318750000e-01", VPFac = "9.69344700036820e-01", KW = "5.12594224560177e-14", KF = "2.46750587115740e-03", KS = "1.07228495518292e-01", K0 = "2.83918818040155e-02", K1 = "1.23204547404511e-06", K2 = "9.14637929020938e-10", FugFac = "9.96810458692103e-01", KB = "2.23064975910959e-09", KP1 = "2.01400858730558e-02", KP2 = "9.31686820250764e-07", KP3 = "1.40307711726248e-09", KSi = "3.62200612635021e-10", KH2S = "2.60000514374855e-07", KNH3 = "4.74226503827862e-10", KCa = "4.92763187414538e-07", KAr = "7.39194729649679e-07")

julia> modern_concs = PALEOcarbchem.calc_modern_default_concs(35.0, Options=options);

julia> map(x -> @sprintf("%.14e", x), modern_concs)
(TF = "6.83258396883673e-05", TS = "2.82354341328601e-02", TB = "4.15700000000000e-04", Ca = "1.02845697008497e-02")

julia> input_concs = (TCi=2000e-6, TS=modern_concs.TS, TF=modern_concs.TF, TSi=1e-3, TP=1e-6,  TB=modern_concs.TB, TH2S=1e-6, TNH3=1e-6, Ca=modern_concs.Ca);  # all in mol kg-1

julia> res = zeros(length(PALEOcarbchem.ResultNames));

julia> (pHfree, steps) = PALEOcarbchem.calculatepHfromTATC!(res, C, options, 2300e-6, input_concs);

julia> @printf("%.14f", pHfree)
8.04695972423890

julia> steps
5

julia> res_NT = NamedTuple{PALEOcarbchem.ResultNames}(res); # convert to NamedTuple for convenience

julia> map(x -> @sprintf("%.14e", x), res_NT)
(pHfree = "8.04695972423890e+00", H = "8.97512024441008e-09", OH = "5.71127974446284e-06", TA = "2.30000000000000e-03", dTAdpH = "6.71952166583155e-04", TS = "2.82354341328601e-02", HSO4 = "2.36333069917159e-09", TF = "6.83258396883673e-05", HF = "2.48522365702671e-10", TCi = "2.00000000000000e-03", CO2 = "1.31351930129921e-05", HCO3 = "1.80311290118315e-03", CO3 = "1.83751905803853e-04", CAlk = "2.17061671279086e-03", fCO2 = "4.62639042514413e-04", pCO2 = "4.64119370418157e-04", xCO2dryinp = "4.78797037215479e-04", TB = "4.15700000000000e-04", BAlk = "8.27503245712347e-05", TP = "1.00000000000000e-06", H3PO4 = "3.68183719676091e-15", H2PO4 = "8.26200822875433e-09", HPO4 = "8.57660283752523e-07", PO4 = "1.34077704336885e-07", PengCorrection = "0.00000000000000e+00", PAlk = "1.12581568874446e-06", TSi = "1.00000000000000e-03", SiAlk = "3.87906357916080e-05", TH2S = "1.00000000000000e-06", H2S = "3.33677816472646e-08", HSAlk = "9.66632218352735e-07", TNH3 = "1.00000000000000e-06", NH4 = "9.49813831954437e-07", NH3Alk = "5.01861680455629e-08", Ca = "1.02845697008497e-02", OmegaCA = "3.83512675291204e+00", OmegaAR = "2.55657840498853e+00")

julia> pHtot = PALEOcarbchem.mappHscale(C, pHfree, Val(3), Val(1), res_NT.TS, res_NT.TF);

julia> @printf("%.14f", pHtot)
7.94544626702045

ComponentsAll::NamedTuple
ComponentsAllStrings::Vector{String}

All available components

julia> PALEOcarbchem.ComponentsAll
(Ci = Val{true}(), S = Val{true}(), F = Val{true}(), B = Val{true}(), P = Val{true}(), Si = Val{true}(), H2S = Val{true}(), NH3 = Val{true}(), Omega = Val{true}())

julia> PALEOcarbchem.ComponentsAllStrings  # all available components
("Ci", "S", "F", "B", "P", "Si", "H2S", "NH3", "Omega")

get_components_inputs(compstoenable) -> (components::NamedTuple, concinputs::Vector{String})

Returns components and concinputs required for specified compstoenable

compstoenable is a list (Vector or Tuple) with component names as Strings, eg ["Ci", "B", "Si"]

Returns:

• a NamedTuple components with those components present in compstoenable as Val{true} and others as Val{false}
• a Vector concinputs of input concentrations required (P requires TP etc).

Examples

All components:

julia> compsall, concinputsall = PALEOcarbchem.get_components_inputs(["Ci", "S", "F", "B", "P", "Si", "H2S", "NH3", "Omega"])
((Ci = Val{true}(), S = Val{true}(), F = Val{true}(), B = Val{true}(), P = Val{true}(), Si = Val{true}(), H2S = Val{true}(), NH3 = Val{true}(), Omega = Val{true}()), ["TCi", "TS", "TF", "TB", "TP", "TSi", "TH2S", "TNH3", "Ca"])

julia> compsall
(Ci = Val{true}(), S = Val{true}(), F = Val{true}(), B = Val{true}(), P = Val{true}(), Si = Val{true}(), H2S = Val{true}(), NH3 = Val{true}(), Omega = Val{true}())

julia> println(concinputsall) # concentrations required (in addition to TAlk)
["TCi", "TS", "TF", "TB", "TP", "TSi", "TH2S", "TNH3", "Ca"]

All optional components disabled (just H2O)

julia> compsminimal, concinputsminimal = PALEOcarbchem.get_components_inputs([])
((Ci = Val{false}(), S = Val{false}(), F = Val{false}(), B = Val{false}(), P = Val{false}(), Si = Val{false}(), H2S = Val{false}(), NH3 = Val{false}(), Omega = Val{false}()), String[])

julia> println(concinputsminimal)
String[]

ConstNames::Tuple

All available CO2SYS effective equilibrium constant names

CNameIdx::NamedTuple

Indices of constants in a Vector of length(ConstNames)

getConstDescription(constName::AbstractString) -> (constName, units, description)

Get human-readable description of constName

calc_constants!(
    Cout, TempC_input, Pdbar, Sal_input; 
    Options=(; WhichKs=Val(10), WhoseKSO4=Val(1), pHScale=Val(3), Components=ComponentsAll))

Calculate carbonate chemistry constants, results are returned in Cout.

Arguments:

• Cout::Vector: (output) Vector of length length(ConstNames) with calculated equilibrium constants
• TempC_input: temperature deg C. Will be limited to valid range for constants in use (from calc_limits)
• Pdbar: pressure, dbar.
• Sal_input: salinity. Will be limited to valid range for constants in use (from calc_limits)
• Options: NamedTuple with fields:
  • WhichKs: Val(Int), default Val(10)
  • WhoseKSO4: Val(Int), default Val(1)
  • pHScale: Val(Int), default Val(3) (free - required for calculateTAfromTCpHfree!)
  • Components: components to enable as generated by get_components_inputs

WhichKs:

K1 K2 dissociation constants that are to be used. Val(Int) where Int is:

• 1 = Roy, 1993 T: 0-45 S: 5-45. Total scale. Artificial seawater.
• 2 = Goyet & Poisson T: -1-40 S: 10-50. Seaw. scale. Artificial seawater.
• 3 = HANSSON refit BY DICKSON AND MILLERO T: 2-35 S: 20-40. Seaw. scale. Artificial seawater.
• 4 = MEHRBACH refit BY DICKSON AND MILLERO T: 2-35 S: 20-40. Seaw. scale. Artificial seawater.
• 5 = HANSSON and MEHRBACH refit BY DICKSON AND MILLERO T: 2-35 S: 20-40. Seaw. scale. Artificial seawater.
• 6 = GEOSECS (i.e., original Mehrbach) T: 2-35 S: 19-43. NBS scale. Real seawater.
• 7 = Peng (i.e., originam Mehrbach but without XXX) T: 2-35 S: 19-43. NBS scale. Real seawater.
• 8 = Millero, 1979, FOR PURE WATER ONLY (i.e., Sal=0) T: 0-50 S: 0.
• 9 = Cai and Wang, 1998 T: 2-35 S: 0-49. NBS scale. Real and artificial seawater.
• 10 = Lueker et al, 2000 T: 2-35 S: 19-43. Total scale. Real seawater.
• 11 = Mojica Prieto and Millero, 2002. T: 0-45 S: 5-42. Seaw. scale. Real seawater
• 12 = Millero et al, 2002 T: -1.6-35 S: 34-37. Seaw. scale. Field measurements.
• 13 = Millero et al, 2006 T: 0-50 S: 1-50. Seaw. scale. Real seawater.
• 14 = Millero et al, 2010 T: 0-50 S: 1-50. Seaw. scale. Real seawater.

WhichKs=Val(10) (Lueker etal (2000) parameters.)

Recommended as 'best practice' by Dickson (2007) (as cited in Orr, J. C., Epitalon, J-M., & Gattuso, J-P. (2014). Comparison of seven packages that compute ocean carbonate chemistry. Biogeosciences Discussions, 11(4), 5327�5397. doi:10.5194/bgd-11-5327-2014 )

WhoseKSO4:

KSO4 dissociation constants that are to be used, in combination with the formulation of the borate-to-salinity ratio to be used. Having both these choices in a single argument is somewhat awkward, but it maintains syntax compatibility with the previous version. Val(Int) where Int is:

• 1 = KSO4 of Dickson & TB of Uppstrom 1979 (PREFERRED)
• 2 = KSO4 of Khoo & TB of Uppstrom 1979
• 3 = KSO4 of Dickson & TB of Lee 2010
• 4 = KSO4 of Khoo & TB of Lee 2010

pHScale: Set pH scale.

Val(Int) where Int is:

• 1 = Total scale
• 2 = Seawater scale
• 3 = Free scale [required for PALEOcarbchem solvers]
• 4 = NBS scale

Implementation

Modified CO2SYS Constants, split into functions for maintainability

Comments from the original code:

% SUB Constants, version 04.01, 10-13-97, written by Ernie Lewis.
% Inputs: pHScale%, WhichKs%, WhoseKSO4%, Sali, TempCi, Pdbar
% Outputs: K0, K(), T(), fH, FugFac, VPFac
% This finds the Constants of the CO2 system in seawater or freshwater,
% corrects them for pressure, and reports them on the chosen pH scale.
% The process is as follows: the Constants (except KS, KF which stay on the
% free scale - these are only corrected for pressure) are
%       1) evaluated as they are given in the literature
%       2) converted to the SWS scale in mol/kg-SW or to the NBS scale
%       3) corrected for pressure
%       4) converted to the SWS pH scale in mol/kg-SW
%       5) converted to the chosen pH scale
%
%       PROGRAMMER'S NOTE: all logs are log base e
%       PROGRAMMER'S NOTE: all Constants are converted to the pH scale
%               pHScale% (the chosen one) in units of mol/kg-SW
%               except KS and KF are on the free scale
%               and KW is in units of (mol/kg-SW)^2

Julia-specific details:

WhichKs, WhoseKSO4, pHScale are passed as Types (using Val(Int)), not integer values, so that Julia can work out which constants and functions to call at compile time. Similarly, Comps is a Type which encodes the selection of components, so the appropriate code is generated (components included or excluded) at compile time.

calc_limits( Val{WhichKs::Int}) -> (TminC, TmaxC, salmin, salmax)

Return range limits for constant set WhichKs:

• 1 = Roy, 1993 T: 0-45 S: 5-45. Total scale. Artificial seawater.
• 2 = Goyet & Poisson T: -1-40 S: 10-50. Seaw. scale. Artificial seawater.
• 3 = HANSSON refit BY DICKSON AND MILLERO T: 2-35 S: 20-40. Seaw. scale. Artificial seawater.
• 4 = MEHRBACH refit BY DICKSON AND MILLERO T: 2-35 S: 20-40. Seaw. scale. Artificial seawater.
• 5 = HANSSON and MEHRBACH refit BY DICKSON AND MILLERO T: 2-35 S: 20-40. Seaw. scale. Artificial seawater.
• 6 = GEOSECS (i.e., original Mehrbach) T: 2-35 S: 19-43. NBS scale. Real seawater.
• 7 = Peng (i.e., originam Mehrbach but without XXX) T: 2-35 S: 19-43. NBS scale. Real seawater.
• 8 = Millero, 1979, FOR PURE WATER ONLY (i.e., Sal=0) T: 0-50 S: 0.
• 9 = Cai and Wang, 1998 T: 2-35 S: 0-49. NBS scale. Real and artificial seawater.
• 10 = Lueker et al, 2000 T: 2-35 S: 19-43. Total scale. Real seawater.
• 11 = Mojica Prieto and Millero, 2002. T: 0-45 S: 5-42. Seaw. scale. Real seawater
• 12 = Millero et al, 2002 T: -1.6-35 S: 34-37. Seaw. scale. Field measurements.
• 13 = Millero et al, 2006 T: 0-50 S: 1-50. Seaw. scale. Real seawater.
• 14 = Millero et al, 2010 T: 0-50 S: 1-50. Seaw. scale. Real seawater.

calc_modern_default_concs(Sal_input, Options=(; WhichKs=Val(10), WhoseKSO4=Val(1)) -> (; TF, TS, TB, Ca)

Calculate modern seawater default concentrations as a function of salinity Sal_input

See calc_constants! for documentation for Options.

Returns

- TF: (mol kg-1) Total Fluorine modern value from salinity
- TS: (mol kg-1) Total Sulphate modern value from salinity
- TB: (mol kg-1) Total Boron modern value from salinity
- Ca: (mol kg-1) Ca modern value from salinity

Example

julia> options = (; WhichKs=Val(10), WhoseKSO4=Val(1));

julia> modern_concs = PALEOcarbchem.calc_modern_default_concs(35.0, Options=options);

julia> map(x -> @sprintf("%.14e", x), modern_concs)
(TF = "6.83258396883673e-05", TS = "2.82354341328601e-02", TB = "4.15700000000000e-04", Ca = "1.02845697008497e-02")

ResultNames::Tuple

All available CO2SYS result names

RNameIdx::NamedTuple

Indices of results in a Vector of length(ResultNames)

getResultDescription(resultName::AbstractString) -> (resultName, units, description)

Get human-readable description of resultName

calculateTAfromTCpHfree!(res, C, Options, concs, pHfree; do_dTAdpH=Val(false)) -> (TA, dTAdpH)

Calculate TAlk, and speciation, given pH and conserved concentrations (total DIC, Si, P, SO4, B, F, ...)

Returns:

• TA: Total Alk, mol/kg-sw
• dTAdpH: derivative, or NaN if do_dTAdpH=Val(false)

Arguments:

• res: (output) Vector res of length length(ResultNames) with details of TA contributions etc
• C: constants from calc_constants!. NB: must be on Free pH scale.
• concs::NamedTuple: (mol kg-sw) total concentrations for sulphate, fluoride, and each optional component of alkalinity enabled in C
• pHfree: pH on free scale

Implementation

modified from Matlab CO2SYS CalculateTAfromTCpH(TCi, pHi) to use Free pH scale

calculatepHfromTATC!(
    res, C, Options, conc_TAx, concs; 
    pHstart=8.0, pHTol=100*eps()
) -> (pHfree, steps)

Call calculateTAfromTCpHfree! to iteratively solve for pH given (i) a starting value and pH tolerance, and (ii) Alk and conserved concentrations (total DIC, Si, P, SO4, B, F, ...)

Intended for use eg in an ocean model to enable a single Newton-Raphson step each model timestep (ie pHstart from previous value in that grid cell, pHtol set to some large number)

Returns:

• pHfree: pH on free scale
• steps: number of Newton iterations

Arguments:

• res: (output) Vector res of length length(ResultNames) with details of TA contributions etc
• C: constants from calc_constants!. NB: must be on Free pH scale.
• conc_TAx: total Alk, mol kg-sw
• concs: other input total concentrations, see calculateTAfromTCpHfree!
• pHstart: starting value (free pH scale)
• pHtol: tolerance (accuracy required)

Implementation

Modified from CO2SYS SUB CalculatepHfromTATC, version 04.01, 10-13-96, written by Ernie Lewis. NB: recoded here to use calculateTAfromTCpHfree!, hence works on free pH scale.

Units: mol / kg-sw

calculateOmega(C, CO3,Ca) -> (OmegaCA, OmegaAR)

Calculate carbonate saturation.

Returns omega, the solubility ratio, for calcite and aragonite. This is defined by: Omega = [CO3–]*[Ca++]./Ksp, where Ksp is the solubility product (either KCa or KAr).

Arguments:

• C - constants from calc_constants!
• CO3 - carbonate ion concentration, mol/kg-sw
• Ca - calcium concentration, mol/kg-sw

mappHscale(C, pHin, scalein, scaleout, TS, TF) -> pHout

Map pH scale at pressure, temperature, salinity defined by constants C

NB: Total, SW scale are not well defined unless using default contemporary values for TS (sulphate) and TF (fluorine)

Arguments:

• C: constants Vector from calc_constants!
• pHin: input pH on scale scalein
• scalein, scaleout :
  • Val(1) = Total scale
  • Val(2) = Seawater scale
  • Val(3) = Free scale
  • not implemented: Val(4) = NBS scale
• TS: total sulphate (mol kg-1)
• TF: total fluoride (mol kg-1)